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Abstract 

We study the dynamics of a reaction-diffusion system composed of two mutually 
coupled excitable fibers. We focus on the situation in which dynamical properties 
of the two fibers are not identical because of the parameter difference between the 
fibers. Using the spatially one-dimensional FitzHugh-Nagumo equations as a model 
of a single excitable fiber, we show that the system exhibits a rich variety of dynam- 
ical behavior, including soliton-like collision between two pulses, recombination of 
a solitary pulse and synchronized pulses, and overtaking of a slow-moving solitary 
pulse by fast-moving synchronized pulses. 
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1 Introduction 



Excitability is a ubiquitous dynamical property encountered in many fields 
of science and plays important roles in the functional aspects of many living 
systems such as transmission of electronic signals in neural and cardiac sys- 
tems [1,2,3]. It has been found that spatially extended excitable media, which 
are modeled in the framework of reaction-diffusion systems, show a rich vari- 
ety of dynamical behavior including propagating pulses and target waves [4], 
spiral waves [5], and spatio-temporal chaos [6]. 
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Problems resulting from intra-medium interactions of spatially localized pat- 
terns in an excitable medium have attracted great interest. One of the distin- 
guishing features of such an interaction is that two propagating pulses in an 
excitable medium annihilate each other upon head-on collision [7]. In general, 
dissipative systems feature not only annihilation of excited waves; they have a 
variety of interactions among spatially localized patterns that behave like an 
elastic object upon collision and scatter in various ways [8,9,10,11,12,13,14,15]. 

On the other hand, inter-media interactions, i.e., interactions among spatially 
localized patterns in multilayered excitable media, should be also of great im- 
portance from the practical viewpoint. For example, in several nerve systems 
such as the hippocampus, olfactory nerves, corpus callosum, spinal column, 
peripheral nerves, and cerebellum, it is observed that huge nerve axons are 
arranged in densely packed bundles so that neighboring neurons can electron- 
ically communicate with each other. Beginning with the pioneering works of 
Katz and Schmitt [16] and Arvanitaki [17], studies on electrical axo-axonal 
interactions have had a long history [18,19,20,21,22,23,24]. Besides nerve sys- 
tems, similar bundled structures are observed in cardiac systems such as the 
bundle of His and Purkinje fibers in the myocardium. Mathematical models 
composed of coupled reaction-diffusion systems [25,26,27,28,29,30,31,32,33,34,35,36,37] 
have been used for elucidating the pattern dynamics in these parallel fibers. 
Coupled react ion- diffusion systems are also used as models of information pro- 
cessing between neural assemblies [38,39]. In non-biological experiments, it has 
been reported that mutual synchronization between two chemical waves occurs 
in the Belousov-Zhabotinsky chemical reaction system with cross-membrane 
coupling [40] and camera- video projector coupling [41]. Understanding and 
controlling pattern dynamics in coupled reaction-diffusion systems are impor- 
tant research subjects in many applications. 

In this paper, we report on the pulse dynamics that emerge from a system 
of two mutually coupled excitable fibers when the dynamical properties of 
two excitable fibers are not identical. Such a situation is not uncommon. For 
example, the diameters of real neuronal fibers are generally not equal. This 
situation is modeled as a difference of diffusion coefficients for each fiber in 
a reaction-diffusion system. Using the spatially one-dimensional FitzHugh- 
Nagumo equations [42,43] as a model of a single excitable fiber, we show that 
in some cases, two propagating pulses do not annihilate upon head-on colli- 
sion and are reconstructed with unchanged spatial profiles like solitons [44]. 
Other interesting and somewhat unexpected pulse dynamics, including recom- 
bination of a solitary pulse and synchronized pulses and overtaking of a slow- 
moving solitary pulse by fast-moving synchronized pulses, are also shown. To 
our knowledge, these pulse dynamics have not yet been reported in previous 
studies on coupled reaction-diffusion systems. 

This paper is organized as follows. A brief description of the mathematical 
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model for coupled excitable fibers is given in Sec. 2. Section 3 introduces 
the reentrant wave that is well-known pattern dynamics observed in coupled 
reaction-diffusion systems, as well as the soliton-like pulse collision that we 
first report in the present paper. Section 4 discusses the stability of the syn- 
chronized pulses as the difference between the intra-diffusion coefficients of 
the two fibers changes. Section 5 provides several examples of pulse dynam- 
ics associated with the destruction of the synchronized pulses. Section 6 is 
a summary and discusses the importance of our findings in the context of 
neuroscience. 



2 Model 



We consider two mutually coupled one- dimensional FitzHugh-Nagumo (FHN) 
fibers. The system consists of the following equations: 



{iii = — a)(l — ui) — v\ + KiV 2 ui + e(u 2 — ui) 

ui = r(tti - 7V1), 

r (1) 

I u 2 — u 2 (u 2 - a)(l - u 2 ) -v 2 + k 2 V 2 u 2 + e(ui - u 2 ) 
{ v 2 = t(u 2 - 7^2). 

Subscripts "1" and "2" denote the first and the second fibers. The state vari- 
ables 1*1,2 = Ui j2 (x,t) and vi j2 = vi i2 (x,t), where x G [0, L] and t G [0, oo) are 
space and time coordinates, are the activators (membrane potentials) and the 
inhibitors (recovery variables), respectively. The parameters of the reaction ki- 
netics are fixed as a = 10" 1 , r = 2 x 10 -3 , and 7 = 2.5 so that a local kinetics 
shows an excitable property; i.e., a small but finite perturbation to the resting 
state (u, v) = (0,0) leads to a large excursion. The terms V 2 Ui t2 = d 2 Ui j2 /dx 2 
represent intra-fiber diffusions and K\ and k 2 are their coefficients. The value 
of K\ is fixed at 0.25 throughout this paper. The mutual interaction between 
two excitable fibers is introduced as linear coupling terms t{u\,2 ~ u 2,i) fo r 
activators. We take e and n 2 to be the control parameters. 

The numerical simulations use the Euler integration scheme with a time step 
At = fO -2 . Diffusion terms at a spatial point Xj(= iAx) with a spatial step 
Ax = 5 x 10 _1 are approximated as V 2 Ui j2 (xi) = (l/(Ax) 2 )(Mi i2 (xj_i) — 
2-u lj2 (xj) + Ui t2 (x i+ i)). Periodic boundary conditions are employed for both 
fibers: ui )2 (0, £) = u 1:2 {L,i) and 1*1,2(0, £) = Vi t2 (L,t). 
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3 Reentrant Wave and Soliton-like Pulse Collision 



Here, we investigate the pulse dynamics of the system of Eqs. (1) when a 
pulse propagating to the right is initiated on fiber 1, and fiber 2 is in the 
global resting state. To prepare these states, we consider the following initial 
conditions: 
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for fiber 1, and 



u 2 {x, 0) = v 2 {x, 0) = Vsce[0,L], (4) 

for fiber 2. Applying these initial conditions and setting the interaction be- 
tween fibers e to zero for t < t , where t is a short interval, we obtain the 
stationary state mentioned above. We redefine these required states as initial 
conditions, and interaction between fibers is taken into account. When Ki and 
k 2 are identical, we observe the following four different phases after the ini- 
tial transient dies out with the increase of inter-fiber coupling strength e: (i) 
a solitary pulse propagating in fiber 1, (ii) formation of the reentrant wave, 

(iii) the global resting state after a finite repetition of reentrant waves, and 

(iv) synchronized pulses propagating in both fibers. In addition to the these 
phases, when Ki and k 2 are not equal, a soliton-like collision occurs depending 
on the control parameters. The following sections deal with these observations 
in detail. 



3.1 Identical Case 



First, we show what happens when K\ and k 2 are identical; i.e., K\ = k 2 = 0.25. 

As shown in Fig. 1, for a sufficiently small value of e, a propagating pulse in 
fiber 1 does not significantly affect fiber 2, and only a sub-threshold excitation 
typically appears as a small amplitude pulse in fiber 2. 
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Fig. 1. Profile of a solitary pulse in fiber 1 with e = 5 x 10~ 3 . A tiny sub-threshold 
excitation is also induced in fiber 2 through the inter-fiber interaction. Arrows in- 
dicate the direction of propagation. 

With increasing e, the qualitative features of the dynamics change dramat- 
ically. When e c i ~ 7.2058 x 10~ 3 , a propagating pulse in fiber 1 induces a 
supra-threshold excitation in fiber 2. Figure 2 shows a series of spatial profiles 
with e = 7.3 x 10~ 3 . The new excitation in fiber 2 splits into two pulses moving 
leftward and rightward as shown in Fig. 2 (a). The right-propagating pulse in 
fiber 2 immediately becomes synchronized with the one in fiber 1, whereas 
the left-propagating pulse in fiber 2 remains solitary (Fig. 2 (b)). Because 
there is a refractory region behind the right-propagating pulse in fiber 1, a 
time interval is needed to develop the subsequent excitation in fiber 1 induced 
by the left-propagating pulse in fiber 2 (Fig. 2 (c)). Consequently, two pulses 
propagating in opposite directions emerge from this subsequent excitation in 
fiber 1, and this new left-propagating pulse is synchronized with the previ- 
ously generated pulse in fiber 2, while the right-propagating one is alone for a 
time. This right-propagating pulse will eventually cause another excitation in 
fiber 2 as shown in Fig. 2 (d). These processes repetitively occur in a specific 
region. The dynamical pattern associated with such repetitions is called the 
reentrant wave [29,30,31,32]. The reentrant wave is considered to be an origin 
of fatal heart diseases such as tachycardias and fibrillation [45,46]. 

In the reentrant wave phase, alternate generation of pulses repeats with a 
characteristic period T. This characteristic period is depicted in Fig. 3 as a 
function of the coupling strength. The period shows a power law with exponent 
1/2, i.e., T ~ |e — edl^ 1 / 2 near the transition point e c \. This result implies that 
the saddle-node bifurcation is the onset mechanism of the reentrant wave. 

The reentrant wave has two qualitatively different features depending on the 
inter-fiber coupling strength e. 1) For smaller values of e, the core position of 
the reentrant wave is fixed in time as shown in Fig. 4 (a). 2) As the coupling 
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Fig. 2. A series of snapshots (a) to (d) for the reentrant wave with e = 7.3 x 10 -3 . 
t = 180, 450, 580, and 960. (a) A supra-threshold excitation is induced in fiber 2 
by a right-propagating pulse in fiber 1, and it splits into two pulses propagating 
in opposite directions, (b) The right-propagating pulse in fiber 2 is synchronized 
with that in fiber 1, whereas the left-propagating pulse remains alone, (c) The 
left-propagating pulse in fiber 2 induces a new supra-threshold excitation in fiber 
1, and it splits into two pulses propagating in opposite directions, (d) In the same 
way, the right-propagating pulse in fiber 1 induces a supra-threshold excitation in 
fiber 2. This alternative excitation in the two fibers is the source of the reentrant 
wave. 

strength becomes larger, the core position of the reentrant wave starts to drift 
to the right. A spatio-temporal plot for the drifted reentrant wave is shown in 
Fig. 4 (b). The direction of drift depends on the initial conditions, and is the 
same as that of the propagating pulse in fiber 1. 

For e > e C 2 ~ 2.83 x 10 -2 , the reentrant wave disappears after a finite repetition 
as shown in Fig. 4 (c). The lifetime of the transient reentrant wave decreases 
as the value of e — e c2 increases. When a solitary pulse propagates in a fiber, 
there is a zone in the other fiber that is stimulated by interacting with the 
pulse. This excitatory induction is stronger for larger e. Thus, the time interval 
required to generate the new pulse becomes shorter, and the period of alternate 
pulse generation decreases. Furthermore, for a given value of e, the period is 
initially longer, and converges to a stationary value that is smaller than the 
initial period. In other words, the distance between the points at which new 
excitations emerge decreases over time. For e C 2 < e < e c3 , a new excitation 
at one time does not split into two propagating pulses because the half side 
of a new excitation meets the refractory region. The repetition subsequently 
terminates. 

Finally, for sufficiently large values of e > e c3 ~ 7.2 x 1CT 2 , the initial pulse 
in fiber 1 generates a new excitation in fiber 2. However, in this case, the 
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Fig. 3. Characteristic period T of the reentrant wave as a function of inter-fiber 
coupling strength e. Inset: The same plot on a logarithmic scale. L=10 3 . The broken 
line in the inset indicates \e — e c i| -1 / 2 . 

new excitation does not form two pulses propagating in opposite directions. 
Only the right-propagating pulse emerges from the excitation and immediately 
synchronizes with the initially generated pulse in fiber l(Fig. 4 (d)). 

To characterize these phases, we use the following spatially coarse grained 
quantity: 



^(t) = \j j I dx(u li2 (x,ty + Vl 4x,ty). (5) 

Figure 5 shows the time averages (cr lj2 (£))t as functions of the inter-fiber cou- 
pling strength e. The four phases of spatio-temporal patterns - (i) solitary 
pulse, (ii) reentrant wave, (iii) transient reentrant wave, and (iv) synchronized 
pulses - are clearly distinguished as {(y\^if))t changes. 

3.2 Nonidentical Case 

Next, turn our attention to the case in which Ki and k 2 are not identical. This 
subsection includes new findings which have not been reported in previous 
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Fig. 4. Spatio-temporal plots for (a) the reentrant wave: e = 1.0 x 10~ 2 , (b) the reen- 
trant wave with drift: e = 2.8 x 10~ 2 , (c) the transient reentrant wave: e = 3 x 10 -2 , 
and (d) synchronized pulses: e = 10 _1 . The red lines indicate the pulses locations 
on fiber 1, which is estimated with the criterion that u\ is above 0.7. The green 
lines indicate the pulse locations on fiber 2. L = 250. 



studies on coupled reaction-diffusion media. In the following, we set Ki to 0.25 
and k 2 to 0.09, and focus on the pulse dynamics in fiber 2. 

As described in the identical case, when e is too small, a sub-threshold exci- 
tation in fiber 2 is generated under the influence of the pulse propagating in 
fiber 1. The solitary pulse does not significantly affect the other fiber. 
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Fig. 5. Time averaged norms (<7i t 2(t))t y s- e for the case of k\ = k%. Stationary states 
are (i) a solitary pulse in the fiber 1, (ii) the reentrant wave, (iii) the global resting 
state after the transient reentrant waves, and (iv) synchronized pulses. L = 10 3 . 

However, when e exceeds e c /i ~ 7.05 x 10 -3 , an excitation in fiber 2 is induced 
by the propagating pulse in fiber 1. This excitation generates two pulses prop- 
agating in opposite directions in fiber 2. This situation occurs in the same 
manner as described in the identical fibers case. For non-identical fibers, asym- 
metrical excitations, which we call one-way excitations, appear in the following 
way. Through inter-fiber interactions, the left-propagating pulse in fiber 2 ac- 
tivates a localized zone in fiber 1, which corresponds to the action potential 
region in fiber 2. Because diffusion generally acts as a smoother for a given 
spatial inhomogeneity, the activated zone in fiber 1 can not be excited due to 
the larger diffusion coefficient (k\ = 0.25). On the other hand, the activated 
zone in fiber 2 is easily excited to form a pulse because of the smaller coefficient 
(k 2 = 0.09), resulting in the one-way excitation. Indeed, compared with the 
case of identical fibers, a stronger infra-fiber coupling is needed to form the 
reentrant wave. The one-way excitation is one of the new phases that emerge 
through the interaction between non-identical fibers, and it exists between 
the solitary pulse and the reentrant wave phases. Furthermore, we found the 
following interesting pulse dynamics in this phase. 



Figure 6 shows typical snapshots of spatial profiles in the one-way excitation 
phase. Due to the periodic boundary condition, two pairs of pulses propagat- 
ing in opposite directions become close to each other (Fig. 6 (a)); i.e., one pair 
is composed of synchronized pulses, which are two supra-threshold excitations 
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Fig. 6. A series of snapshots (a) to (d) for the soliton-like pulse collision with 
«i = 0.25, k 2 = 0.09, and e = 8 x 10" 3 . t = 1.0 X 10 3 ,1.3 x 10 3 ,1.5 x 10 3 , and 
1.9 x 10 3 . (a) Head-on collisions occur in both fibers, (b) Mutual annihilation oc- 
curs in fiber 2 whereas the pulse prevails over the sub-threshold excitation in fiber 1. 
(c) A supra-threshold excitation that splits into two propagating pulses is induced 
in fiber 2 by the pulse in fiber 1. (d) The spatial profiles of all pulses recover after 
head-on collisions. 



both in fibers 1 and 2 with a tiny delay, and the other pair is composed of 
a solitary pulse propagating in fiber 2 and a sub-theshold excitation in fiber 
1. Eventually, head-on collisions occur in both fibers (Fig. 6 (b)). In fiber 1, 
the head-on collision between the supra-threshold pulse and the sub-threshold 
excitation does not lead to mutual annihilation; i.e., the sub-threshold excita- 
tion vanishes and the right-propagating pulse persists. On the other hand, in 
fiber 2, the two propagating pulses mutually annihilate once by the head-on 
collision. A new supra-threshold excitation, however, is regenerated through 
induction from the right-propagating pulse in fiber 1. This excitation splits into 
two pulses propagating in opposite directions (Fig. 6 (c)). Finally, the spatial 
profiles of all pulses completely recover after the head-on collision (Fig. 6 (d)). 
These kinematic features are similar to the characteristics of solitons [44] in 
integrable systems. A spatio-temporal plot of this solition-like collision is also 
shown in Fig. 7. 

Theoretical studies have predicted that head-on collisions between propagating 
pulses do not always lead to mutual annihilation, but form a variety of nontriv- 
ial dynamical behavior, including soliton-like collisions [8,9,10,11,12,13,14,15], 
and there are experiments that verify their existence [47,48]. 

It is also known that a transition from annihilation to crossing of the pulses 
upon a head-on collision occurs in an FHN system [49] that considers diffusion 
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Fig. 7. Spatio-temporal plot for a soliton-like pulse collision with 
k x = 0.25, k 2 = 0.09, and e = 8 x 10~ 3 . L = 250. 

of the inhibitor v and bistable reaction kinetics; i.e., coexistence of a stable 
limit cycle and a stable fixed point. In real nerve systems, however, diffusion 
for the inhibitor generally does not exist, because the inhibitor represents 
potential-dependent gating variables, and a fiber usually has a simple excitable 
property. In our system, we do not require both diffusion for the inhibitor and 
multi-stable kinetics. Instead, fiber 1 plays the role of an auxiliary field for the 
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Fig. 8. Time averaged norms (oi^i))* vs. e for the case of K\ ^ ki- Inset: An 
enlargement of the region e G [5 x 10 _3 ,10 -2 ]. Six different phases exist in the 
case of coupled non-identical fibers: (i) a solitary pulse in fiber 1, (ii) a soliton-like 
pulse collision, (iii) the reentrant wave, (iv) synchronized pulses after the transient 
reentrant waves, (v) the global resting state after the transient reentrant waves, and 
(vi) synchronized pulses. K\ = 0.25, &2 = 0.09 (Ak = 0.16), and L = 10 3 . 



occurrence of the soliton-like behavior in fiber 2. 



As e increases, the reentrant wave phase and the synchronized pulses phase 
appear in non-identical fibers, as well as in the identical fibers. The phase 
diagram of the nonidentical case is depicted in Fig. 8. For e > e c / 2 , i.e., after 
the formation of reentrant wave, the qualitative behaviors are similar to those 
of the indentical fibers as shown in Fig. 5. However, a new phase appears, as 
indicated by (iv) in Fig. 8. In this phase, synchronized pulses form after a finite 
repetition of the reentrant waves. We also investigated the phase diagram in 
two-dimensional parameter space (e, An), where An = K\ — k^. The result 
is shown in Figure 9 . There are clear bifurcation curves that discriminate 
different phases. The soliton-like phase still survives for tiny Ak, and the 
interval of e for it extends as Ak increases. 
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Fig. 9. Phase diagram in the plane of (e, Ak), Ak = ki — K2- «i is set to 0.25 and 
L = 10 3 . 



4 Stability of Synchronized Pulses 



In the previous section, we considered the case where a right-propagating pulse 
is initiated only in fiber 1 at t = 0. As shown in Figs. 5 and 8, for such initial 
conditions, synchronized pulses are observed for e > e C 3 in the identical case 
and [e G (e c / 3 , e^)] V [e > e c / 5 ] in the nonidentical case. It is expected, however, 
that when Ak = K\ — k 2 is small, synchronized pulses will be in a stable state 
for small e if a pair of synchronized pulses are taken as initial conditions. 
From the viewpoint of information processing in the brain, the stability of 
synchronized propagating pulses in a nerve bundle is of interest. Thus, in this 
section, we consider the stability of the synchronized pulses against changes 
in Ak and e when a pair of synchronized pulses are initiated at t — 0. K\ is set 
to 0.25 and k 2 is varied within an interval [0, 0.25] in the following. To detect 
whether the synchronization is stable or not, we consider the velocities c\ and 
C2 of the two pulses propagating in fiber 1 and 2, respectively. 

Figure 10 (a) plots c\ and c 2 as functions of Ak for a relatively larger e(= 
8 x 10~ 3 ). Although both velocities decrease with increasing Ak, the result of 
c\ = c 2 for any Ak indicates that two propagating pulses remain synchronized. 
On the other hand, for smaller e(= 6 x 10~ 3 ), there is a sudden dip in the graph 
of c 2 vs. Ak at Ak c ~ 0.1 whereas the graph of c\ vs. Ak is continuous, as 
shown in Fig. 10 (b). This result indicates that the synchronization becomes 
unstable for Ak > Ak c . Indeed, even if synchronized pulses are initially given, 
the pulses become de-synchronized and propagate with different velocities. 

The dependence of the critical difference Ak c on e is shown in Figure 11. There 
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Fig. 10. Velocities c\ and C2 of propagating pulses in fibers 1 and 2 as functions of 
An = k% — K2 with (a) e = 8 x 10~ 3 , and (b) e = 6 x 10~ 3 . Open circles indicate 
ci and crosses c%. K\ is set to 0.25. Velocities of solitary pulses in the two fibers 
without the inter- fiber interaction (i.e., e = 0) are plotted as broken lines in (b). 
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Fig. 11. Critical difference Ak c between intra-diffusion coefficients as a function of e. 
K\ is set to 0.25. Beyond the critical value, the synchronous pulses become unstable, 
and split into two solitary pulses with different propagating speeds. 



is a critical inter-fiber coupling strength ~ 6.8 x 10 3 , beyond which the 
synchronization is sustained for any An < 0.25. 
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5 Recombination and Overtaking 



In the soliton-like phase, two pulses in fiber 2 face each other as shown in 
Fig. 6 (a). Here, we consider an initial state such that there are two pulses 
propagating in the same direction in fiber 2; one is synchronized with a pulse 
in fiber 1, and the other is a solitary pulse as depicted in Fig. 12 (a). This 
initial state is also established by an appropriate initial stimulation at the 
same parameter values for which the soliton-like collision is observed. Here, 
we set K\ = 0.25, k 2 = 0.09, and e = 8 x 10 -3 . Starting from the initial state 
(Fig. 12 (a)), we observe interesting dynamical behavior associated with the 
loss of synchronized pulses. 

In the following, we use the notation PI, P2A, and P2B to indicate the three 
excited pulses in Fig. 12 (a). All three pulses propagate to the right, and PI 
and P2A are synchronized with each other. As shown in Fig. 10 (b), the syn- 
chronized pulses are faster than the solitary pulse in fiber 2. Therefore, P2A 
is faster than P2B, and the distance between P2A and P2B decreases with 
time. However, when P2A comes sufficiently close to P2B, the highly concen- 
trated region of the inhibitor behind P2B decelerates P2A, leading to loss of 
synchronization between PI and P2A (Fig. 12 (b)). Since PI is faster than 
P2A and P2B, synchronization between PI and P2B is eventually achieved 
(Fig. 12 (c)). These new synchronized pulses move away from P2A (Fig. 12 
(d)). We call this series of dynamic processes recombination of a solitary pulse 
and synchronized pulses. A corresponding spatio-temporal plot is shown in 
Fig. 13. 

To clarify how this recombination depends on An, we investigate the time 
evolution of the distance I between P2A and P2B in fiber 2 for e = 7.3 x 10 -3 . 
Depending on the value of Ak, two qualitatively different results are obtained, 
as shown in Figs. 14 (a) and (b). 

Figure 14 (a) plots I vs. t for four different values of Ak. For each graph, the 
distance between P2A and P2B becomes short over time and approaches a 
local minimum. The time at which I reaches its local minimum corresponds 
to the moment when synchronization between PI and P2A is broken. After 
that moment, there is a time interval in which I slowly increases. This time 
interval corresponds to the recombination process. Focusing on the change in 
the graphs from (i) to (iv) in Fig. 14 (a), it is seen that the elapsed time re- 
quired for recombination gradually increases with decreasing Ak, but diverges 
at Ak = 8.5 x 10 -2 . This result implies that the distance between P2A and 
P2B converges to a stationary value and the recombination is abandoned along 
the way as shown in Fig. 15. This profile is similar to the one in Fig. 12 (b); 
however, all three pulses propagate with the same speed, thus, the relative 
positions among three pulses are fixed in time. We call these dynamics locking 
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Fig. 12. A series of snapshots (a) to (d) for recombination of synchronized and 
solitary pulses with k\ = 0.25,^2 = 0.09, and e = 8 x 1CT 3 . L = 700. 
t = 4.8 x 10 3 ,5.7 x 10 3 ,6.0 x 10 3 , and 6.6 x 10 3 . (a) Synchronized pulses com- 
posed of PI and P2A comes close to P2B. (b) Synchronization is broken by the 
inhibitor behind P2B. (c) New synchronized pulses composed of PI and P2B form, 
(d) After recombination, the new synchronized pulses move away from P2A. 

of propagating pulses. 

For larger values of A/c, the qualitative feature of the graph, I versus t, changes 
as shown in Fig. 14 (b); i.e., the time interval required for the recombination 
becomes very short and a clear V-shaped structure appears. Figure 16 shows 
a series of profiles at An = 0.24. It seems that the slow-moving solitary pulse 
is overtaken by the fast-moving synchronized pulses. 



6 Summary and Discussion 



We investigated pulse dynamics in mutually coupled excitable fibers and showed 
that new pulse dynamics occur: solition-like pulse collisions, whereby the spa- 
tial profiles of pulses completely recover after head-on collisions, recombina- 
tion, whereby the combination of synchronized pulses changes into a new one, 
locking, whereby the recombination process is abandoned along the way, and 
overtaking, whereby a slow-moving solitary pulse is overtaken by fast-moving 
synchronized pulses. These exotic behaviors come from the interaction between 
fibers with different intra-diffusion coefficients. In other words, one fiber acts 
as a "hidden variable" for the other fiber. Thus, excitable media, which have 
FHN-type reaction kinetics, show these exotic behaviors. 

Before ending, we should discuss the significance of the present study in the 



16 



[x 30 time units] 




Fig. 13. Spatio-temporal plot for the recombination of synchronized and solitary 
pulses with k\ = 0.25, K2 = 0.09, and e = 8 x 10~ 3 . L = 500. Snapshots are also 
shown in Fig. 12. 




5000 10000 15000 2000 4000 6000 8000 

t t 

Fig. 14. (a) Time evolution of the distance I between two pulses P2A and P2B in 
fiber 2 (see Figs. 12 and 13). (i) Ak = 0.1, (ii) An = 0.09, (iii) Ak = 0.087, and 
(iv) Ak = 0.085. For smaller Ak, the recombination comes to halt, and a locking 
state appears (see Fig. 15 for snapshot), (b) The graph of I vs. t is also shown for 
Ak = 0.24. For larger values of Ak, the qualitative feature changes; i.e., the behavior 
of fast-moving synchronized pulses overtaking a slow-moving solitary pulse emerges 
(see Fig. 16 for the series of snapshots), e = 7.3 x 10~ 3 and L = 700. 

context of neuroscience. As mentioned in Sec. I, there exist densely packed 
bundles of nerve axons in several regions of the brain. In such a structure, 
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Fig. 15. A snapshot for the phase locked pulses. An = 0.085 («i 
K2 = 0.165) and e = 7.3 x 10~ 3 . The pulse formation is fixed in time. 
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Fig. 16. A series of snapshots (a) to (d) for fast-moving synchronized pulses over- 
taking a slow-moving solitary pulse for An = 0.24 («i = 0.25 and k>2 = 0.01). 
e = 7.3 x 10~ 3 , and L = 700. t = 2.4 x 10 3 , 3.0 x 10 3 , 3.2 x 10 3 , and 3.9 x 01 3 . 



one can expect an electrical interaction between adjacent nerve axons. We 
introduced diffusive couplings 6(1*1,2 — ^2,1) in Eqs. (1) as electrical interactions 
between two excitable fibers. 



In general, however, most nerve axons are insulated by a lipid material called 
the myelin sheath with periodic gaps of exposure called the nodes of Ranvier 
in order to increase the speed and the reliability of conduction of pulses [2]. 
Because the myelin sheath reduces current flow that leaks out across the mem- 
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brane, the electrical interaction between neighboring axons surrounded by 
myelin sheaths is very weak; hence it is said that the effect of electric cou- 
pling between axons is negligible [50]. Theoretical studies [36,37] have also 
shown that the relative locations of the nodes of Ranvier on two neighboring 
myelinated nerve axons influence the degree of synchronization between prop- 
agating pulses. It would be interesting to extend our study to incorporate a 
spatially inhomogeneous structure such as the myelin sheath and the nodes of 
Ranvier into the system of Eqs. (1) for actual biological applications. Such a 
property could be an interesting topic of future work. 

On the other hand, it is also known that there are nerve bundles composed 
of unmyelinated axons. In fact, olfactory nerve axons are unmyelinated and 
arranged in densely packed bundles [23] , and the question of whether there are 
gap junctional couplings between adjacent axons has been also discussed [24]. 
We have mainly focused on the pulse dynamics when the diffusion coefficients 
K\ and K2 in Eqs. (1) are not equal. The difference between coefficients might 
correspond to the difference in the diameters of the two nerve axons. Such 
an assumption is quite natural when considering real nerve axons. Thus, real 
neural systems should have the ability to produce a wide variety of pulse 
dynamics similar to the ones shown in the present paper. 

Our results strongly suggest that a nerve axon is not only a cable for trans- 
mitting an input signal to the terminal but also a functional element that can 
process input information by "crosstalk" between adjacent axons in a bun- 
dle [51]. Such crosstalk is not desirable for the purpose of precise conduction 
of the input pulse, and may be also involved in serious neural diseases such as 
Multiple Sclerosis (MS), which is caused by the loss of the myelin sheath [37]. 
It is also known that the electrical interaction between lesioned nerve axons 
is an origin of neuropathic pain such as hyperalgesia and allodynia [52] . 

On the other hand, crosstalk among nerve axons may play a constructive role 
in the brain. In [23] , it is suggested that electrical interactions between neigh- 
boring axons in the olfactory nerves contribute to olfactory discrimination 
by modulating the frequency of the action potentials in neighboring axons 
and by enhancing the degree of synchronization between their firings. Such a 
situation can be also found in mossy fibers in the hippocampus. The struc- 
ture of hippocampal mossy fibers is similar to that of the olfactory nerves; 
i.e., the nerve axons are unmyelinated and arranged in densely packed bun- 
dles. Thus, electrical axo-axonal interactions may occur in mossy fibers and 
are likely to impact the mechanism of neuronal integration [53]. A recent ex- 
perimental study revealed that depolarization at an axon terminal strongly 
influences the efficiency of synaptic transmission in mossy fibers [54]. In our 
simple mathematical model, a pulse in a fiber influences the propagation of 
pulses in the other fibers. Indeed, synchronization and locking between prop- 
agating pulses are established depending on the parameter, i.e., the coupling 
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strength between fibers, and the difference between the infra-diffusion coeffi- 
cients. The interspike interval in the early stage of propagation will change as 
the pulse travel in fibers. Thus, synchronization and locking regulate the inter- 
spike interval. This fact implies that synaptic transmission might be affected 
by electrical axo-axonal interactions in real neuronal systems. 

Various types of pulse dynamics caused by axo-axonal interactions in the 
nerve bundle may act as information processing in the brain [51]. The results 
presented in this paper provide a new insight into the functional aspects of 
nerve axons. 
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